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Abstract. We present a new hybrid lattice-Boltzmann and Langevin molecular 
dynamics scheme for simulating the dynamics of suspensions of spherical colloidal 
particles. The solvent is modeled on the level of the lattice-Boltzmann method while 
the molecular dynamics is done for the solute. The coupling between the two is 
implemented through a frictional force acting both on the solvent and on the solute, 
which depends on the relative velocity. A spherical colloidal particle is represented 
by interaction sites at its surface. We demonstrate that this scheme quantitatively 
reproduces the translational and rotational diffusion of a neutral spherical particle in 
a liquid and show preliminary results for a charged spherical particle. We argue that 
this method is especially advantageous in the case of charged colloids. 
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1. Introduction 

Understanding the dynamics of colloidal dispersions has been a long-standing problem 
in condensed matter physics. Despite the continuous need for accurate theoretical 
predictions concerning particle mobilities or sedimentation rates, the development of the 
field meets a nearly unpassable barrier. The many-body character of the hydrodynamic 
interactions (HI) between the colloidal particles (i. e. their correlated motion as a result 
of solvent-mediated fast momentum transport) poses a serious challenge for analytic 
theory. On the other hand, the rapid growth of available computer power in the last 
decades has made it more and more feasible to descend to a more basic description 
of colloidal systems rather than trying to develop an adequate macroscopic description. 
Such a step was made recently for the statics of charged colloids, where simulations of the 
primitive electrolyte model, taking the particle character of the charges and fluctuations 
in the charge distribution fully into acount, resolved some important strong correlation 
issues and thus have resulted in a significant refinement of the existing theories of 
screening pQ. In a similar fashion, many-body effects in the HI of suspensions of only 
moderate density are quite important. Although their analytical form is, in principle, 
known in terms of a multipole expansion |2j, it is numerically very cumbersome to take 
these higher-order many-body terms into account in a Brownian or Stokesian dynamics 
[3] simulation, where only the positions of the colloidal particles enter. For this reason, 
it is more convenient and for large particle numbers also more efficient, to take these 
effects into account by simulating the solvent degrees of freedom, and, in particular, 
the momentum transport through the solvent explicitly. An additional bonus is that 
retardation effects, if present, are automatically taken into account. 

Standard Molecular Dynamics (MD), where just Newton's equations of motion 
are integrated numerically, is of course able to simulate HI correctly. However, it is 
very time-consuming to simulate the motion of a huge number of particles explicitly in 
such detail. The time step is governed by the local oscillations of the solvent particles 
in their temporary "cages", which is much faster than the motion of the colloidal 
particles. This information is however not needed for correctly reproducing HI. One 
just needs an appropriate mechanism for momentum transfer on somewhat larger time 
scales (however still small compared to the motion of the colloidal particles). This can 
be done in various ways: Dissipative Particle Dynamics (DPD) [4; is essentially MD, 
where however the particles are made very soft in order to increase the time step, and 
where a momentum-conserving Langevin thermostat is added [5 . Another possibility is 
multi-particle collision dynamics 0, where the solvent is modeled as an ideal gas, and 
collisions are modeled as simple stochastic updating rules which conserve energy and 
momentum. Grid-based methods can be either the direct solution of the Navier-Stokes 
equation by a finite difference method, or the Lattice Boltzmann (LB) [Z| approach, 
where essentially a linearized Boltzmann equation is solved in a fully discretized version 
(i. e. space, time, and velocities are discretized). 

In our opinion, all of these approaches have their advantages and disadvantages. 
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Nevertheless, it seems that all of them are roughly equivalent with respect to their 
computational efficiency, and with respect to the level of accuracy of describing the 
physical phenomena. Therefore, we believe that the choice of method is mainly a matter 
of convenience. 

For our model we have chosen the LB route, mainly because this algorithm is 
particularly simple to implement and parallelize. Moreover, we were inspired by previous 
highly successful simulations of hard-sphere colloidal suspensions |Hl HH US HH H2 HS1 
HH [13 HE] • A further advantage of a grid-based method is that thermal fluctuations 
can be both turned on and off, depending on the physical situation under consideration 
(in a particle method they are always present). This approach is essentially a hybrid 
method, where the solvent is modeled via LB, while the colloidal particles are run with 
MD. 

The original approach by Ladd [HI El EH HU H2j models the colloidal particles as 
extended hollow spheres, while stick boundary conditions at the surface are implemented 
(roughly spoken) via bounce-back collision rules. One should note that for moving 
spheres this involves some minor fluctuations in the fluid mass, which is included within 
a sphere: The moving shell incorporates some new fluid at the front, while it releases 
mass at the rear. The fluctuations are a natural effect of the thermal density fluctuations 
of the fluid, and of the lattice discretization. Using this method, the translational and 
the rotational dynamics of colloidal spheres were accurately described. 

This model has been recently extended to the case of charged systems [T7|, where 
the colloidal particles carry a charge, while the counterions and salt ions are taken into 
account as LB populations, such that the electrostatics is essentially taken into account 
on the level of the Poisson-Boltzmann equation. This method has two disadvantages: 
firstly, the discrete nature of the ions and correlations beyond the Poisson-Boltzmann 
level are not taken into account; secondly, one cannot avoid a leakage of charge into 
(and out of) the sphere (as in the case of mass), such that it is hard (if not impossible) 
to maintain a well-defined Debye layer of charges around it ^H]- For these reasons, it 
is advisable to take the counterions explicitly into account. The disadvantage, however, 
is that only rather small size ratios (size of colloidal particle vs. size of ion) are easily 
accessible. Nevertheless, we believe it is useful to study such a system, with respect 
to both its equilibrium and its nonequilibrium properties. The purpose of the present 
paper is to describe first steps in the development and validation of the corresponding 
model. 

For the coupling of the small ions to the LB hydrodynamics, one would like to 
simply model the former as point particles. Fortunately, such a coupling has been 
recently developed by Ahlrichs and Diinweg (THl ED] with the purpose of studying 
the dynamics of polymer solutions [21]. Each Brownian point particle is assigned 
a phenomenological friction coefficient (, and the coupling is implemented just as 
a dissipative force F = — ({V — u) acting on the particle, where V is the particle 
velocity, while u is the solvent flow velocity at the particle's position, obtained via 
linear interpolation from the surrounding lattice sites. Furthermore, the LB variables 
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at these sites are adjusted to ensure momentum balance, and thermal fluctuations are 
added to both types of degrees of freedom. For more details on implementation and 
validation of this method, see Refs. (THl ED] • 

It is then convenient to use this coupling not only for the ions, but also for the 
colloidal particles, whose larger size is taken into account by modeling them as some 
arrangement of interaction sites. For reasons of efficiency, we take only points at the 
surface of the colloid. Furthermore, the coupling constant £ is chosen rather large, 
which approximates stick boundary conditions. It should be noted that the dissipative 
nature of the coupling does not model any squeezing-out of solvent, hence the sphere is 
filled with the same amount of solvent regardless if it is hollow or if it were filled with 
particles. The solvent within the sphere follows the motion of the surrounding shell 
(with some minor time lag, see below), and therefore one can view the solvent inside 
the sphere as just belonging to the colloidal particle. Adding further particles in the 
sphere's volume would have no effect except coupling the fluid within the sphere even 
tighter to its motion. Within our desired level of accuracy, this turned out not to be 
necessary. The long-time motion of the sphere is the same as that of a corresponding 
hard sphere, as we will show in the present paper. 

We analyze basic dynamic properties of such a model colloid, where we restrict 
ourselves to the case of a single sphere. A detailed analysis is done for the neutral 
case, while some preliminary data for the case of a charged sphere are presented. We 
show that the (neutral) model exhibits the essential dynamical features of a spherical 
colloidal particle in a liquid. In our view, our method provides comparable efficiency 
to Ladd's approach, while being quite straightforward to implement. Furthermore, it 
provides substantial flexibility with respect to the properties of the colloidal surface, 
namely, deformable, permeable, and non-stick surfaces can be easily simulated. 

The remainder of this article is organized as follows: In Sec. 121 we describe our 
simulation model, while Sec. El contains the numerical results on translational and 
rotational diffusion. Finally, Sec. H] concludes with a brief summary. 

2. Model 

Our hybrid simulation method involves two subsystems: the solvent that is modeled via 
LB with fluctuating stress tensor (i. e. we run a constant-temperature version of the 
LB method) and a Langevin MD simulation for the particles immersed in the solvent. 
The LB simulation is performed using the 18- velocity model [H], using the protocol 
described in [13120]. The fluid simulation consists of collision and propagation steps, 
the former being performed with inclusion of the momentum transfer from the solute 
particles (surface beads, and, for charged systems, ions). 

The colloidal particle is represented by a two-dimensional tethered bead-spring 
network consisting of 100 beads, which is wrapped around a ball of a radius a cs (for 
notation, see below), so that the whole construction resembles a raspberry (see Fig. [TJ. 
The network connectivity is maintained via finitely extendable nonlinear elastic (FENE) 
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Figure 1. Raspberry- like model of a colloidal sphere. There is a central large bead of 
radius R — 3 and charge Z = 10. The small beads of radius 1 are connected with their 
nearest neighbors on the surface via FENE bonds. A repulsive soft-core potential is 
also operating between all the monomers. The counterions are moving freely in space 
and interact with the central bead via the Coulomb potential and with the surface 
beads via the repulsive LJ potential. IMoviel 



springs, 

lWr) = -^ln j, (1) 

where k is the spring constant, and Rq the maximum bond extension. Furthermore, the 
beads repel each other by a modified Lennard- Jones (LJ) potential 



V LJ (r) = ^ VrJ Vr J 4, "V (2) 

[0 r> 2 l ' & a ij . 

An additional repulsive LJ bead is introduced at the center of the sphere in order to 
maintain its shape. In Eq. i,j denote either a central ("c") or a surface ("s") bead. 
The unit system is completely defined by the surface bead parameters by setting e ss , 
a ss , and the surface bead mass m s to unity. The interaction between the central bead 
and the surface beads is described by a cs = 3, which is thus the sphere radius, and 
e cs = 8. Furthermore, the FENE spring constant for the surface beads is k = 300 and 
the maximum bond extension is Rq = 1.25. To simulate a charged colloidal particle, 
we place the charge at the central bead, and add an appropriate number of counterions 
(LJ beads with "s" properties) outside the sphere. The electrostatic interaction is taken 
into account via the Coulomb potential 

V el (r) = \ B h B T™i (3) 



A new model for simulating colloidal dynamics 



6 



between the various charges, where the standard Ewald summation technique j22] is 
applied. In Eq. El A B = e 2 / (47ie e r kBT) is the Bjerrum length, k B the Boltzmann 
constant, qi the charge of species i in units of the elementary charge e, and T the 
temperature. 

The LB lattice constant is chosen as one (in our LJ unit system), and the fluid is 
simulated in a cubic box with periodic boundary conditions. The force between the LB 
fluid and the surface beads (or ions) is given by 

F=-((y-u)+f. (4) 

Here, ( is the "bare" [20] friction coefficient, V and u are the velocities of the bead and 
the fluid (at the position of the bead), respectively, while / is a Gaussian white noise 
force with zero mean, whose strength is given via the standard fluctuation-dissipation 
theorem [IQ 120] to keep the surface beads and ions at the same temperature as the 
solvent. The central bead is not coupled to the solvent (here £ = / = 0); as discussed 
in the Introduction, the behavior of the model would change only marginally if such 
a coupling were included. In our simulation we used a friction constant ( = 20, a 
temperature k B T = 1, a fluid mass density p = 0.85, and a kinematic viscosity v — 3, 
resulting in a dynamic viscosity n = 2.55. At least 20000 MD steps were performed 
to equilibrate the initial random bead configuration before the interaction with the LB 
solvent was turned on. Further details on the method can be found in TOi 20J. 

3. Results 

We test the simulation method against basic relations for an isolated sphere in solvent. 
First, we look at the center of mass' velocity relaxation. The simplest experiment to 
perform is a "kick". The sphere is placed in a LB fluid at rest (i. e. without thermal 
noise), and at time t = all particles of the sphere are assigned an identical velocity 
V — 1 in x direction. Fig. El monitors the time behavior of the sphere's center of 
mass velocity, normalized by the initial value. According to linear response theory, 
this relaxation function must be identical to the normalized center-of-mass velocity 
autocorrelation function for Brownian motion in thermal equilibrium, if the initial kick 
is weak enough. This is indeed satisfied, as a comparison of the two curves in Fig. |2] 
shows. For the experiment in thermal equilibrium, we performed 10 runs with different 
random number generator initializations, in order to reduce the statistical uncertainty. 

It is well-known that simulations of Brownian motion in a hydrodynamic solvent 
are always strongly affected by finite size effects. The diffusion constant, and therefore 
also the relaxation function, depend on the linear system size L due to hydrodynamic 
interactions with the periodic images. The diffusion constant exhibits a finite-size 
correction of order R/L |20| |2*H]. Asymptotic behavior can therefore only be expected 
for R/L <C 1, and this is why we performed the experiment in a rather large box of size 
L = 80. For the same reason, the equivalence between "kick" experiment and Brownian 
motion will only hold if the comparison is done for the same box sizes. 
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Figure 2. Normalized translational velocity of the center of mass of the colloidal 
sphere in a "kick" experiment at kgT = 0, and its normalized velocity autocorrelation 
function for Brownian motion at ksT = 1. The dotted curve shows the exponential 
decay derived from the Stokes law and the dashed curve the expected long-time 
asymptotic behavior. 




Figure 3. Mean-square displacement of the center of mass of a neutral sphere at 
ksT = 1 and different simulation box sizes. 
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Figure 4. The long-time self-diffusion coefficient and the rotational diffusion 
coefficient of the colloidal sphere (normalized by the asymptotic Stokes-Einstein and 
Stokes-Einstein-Debye values) as functions of the inverse size of the primary simulation 
cell. 

One clearly sees that the initial decay of the relaxation function is characterized by 
two relaxation processes, one initial fast decay followed by a somewhat slower relaxation. 
Qualitatively, this may be explained as follows: A compact sphere of radius R and 
mass M should exhibit a velocity relaxation which, after transient ballistic motion, 
is initially characterized by an exponential decay exp (— t/r), with a relaxation time 
t = M/(tot- Here, (tot is the total friction coefficient, which we estimate via Stokes' 
law for stick boundary conditions as Qot ~ QirnR ~ 144. However, the effective mass 
is time- dependent: While initially only the mass of the beads M = 101 contributes, at 
later times the fluid within the sphere is dragged as well, such that then the mass is 
roughly estimated asM^ 101 + 4irpR 3 /3 ~ 214. This gives rise to the initial and final 
relaxation times rj n ~ 0.7 and r/j n ~ 1.5. However, the initial relaxation time cannot be 
observed, since in the extreme short-time regime ballistic effects play a role. Conversely, 
the decay with Tf in is clearly visible (see Fig. |2j). 

After t ~ 1, the famous long-time tail [2U ESj (normalized relaxation function 
V(t)/V(0) = Bt~ 3 / 2 ) sets in. The physical mechanism of this slow relaxation is the 
fact that the momentum is conserved and hence transported away diffusively from the 
particle [22] • For this reason, it can only be observed in a sufficiently large box. The 
prefactor of the power law is known in the colloidal limit, B = (l/12)(iVm/p)(7rz/)~ 3 / 2 , 
where N is the number of beads, and m the bead mass HZj- It should be noted 
that the prefactor of the power-law decay of the unnormalized velocity autocorrelation 
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Figure 5. Decay of angular velocity of the colloidal sphere. The different curves 
are marked by the primary simulation box sizes. The solid curves show the expected 
exponential decay according to the Debye law and the long-time asymptotic behavior. 



function does not depend on the properties of the sphere at all, but only on the 
temperature, and the hydrodynamic properties of the solvent I2Z|- The mass 
dependence of B results only from the normalization, i. e. the value of (V 2 ) according 
to the equipartition theorem, (V 2 ) = 3/c#T '/(Nm). From this consideration it is clear 
that the short-time value of the sphere's mass (Nm) enters (at t = the shell and the 
inner fluid are not yet coupled). For our simulation parameters, we find B = 0.342. As 
Fig. 121 shows, the data exhibit the expected behavior very nicely. 

Figure El illustrates the finite size effect in the diffusive properties by plotting 
the mean square displacement of the colloidal sphere as a function of time, for two 
different box sizes L = 6.33.R and L = 20R, for Brownian motion at k^T = 1. In the 
long-time regime, the slope is given by 6D , where D s is the self-diffusion coefficient. 
The figure clearly shows that D s increases with box size. This can be explained in 
terms of hydrodynamic interactions with the periodic images, or, equivalently, in terms 
of suppression of long- wavelength hydrodynamic modes [201 ES| We have therefore 
performed a systematic finite-size analysis, and measured D s for various box sizes by 
integrating the velocity autocorrelation function over time [25] • In Fig. Hwe plot D s /D 
(D denoting the asymptotic Stokes-Einstein value D = k B T / (QirnR) = 6.9 x 10~ 3 ) as 
a function of 1/L; the expected linear behavior [22] is clearly seen. 

We now look at the relaxation of the rotational motion. We performed a similar 
kick experiment as described above, where now an initial angular velocity Uq = 1 was 
provided to the sphere. The data for the normalized decay function u(t)/u(0) are 
presented in Fig. El The sphere dynamics shows a characteristic "raw egg" (damped) 
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rotation pattern with a fast initial decay and a subsequent slower one. For the rotational 
relaxation of a hard sphere with moment of inertia I and rotational friction coefficient 
8iir]R 3 , theory predicts a decay according to the Debye law u(t)/ui(Q) = exp(— t/r r ), 
where the relaxation time is given by r r = I/(8irr]R 3 ) [2H]- Similarly to the effective 
mass, the effective moment of inertia is expected to increase as a function of time, due 
to the time-delayed dragging of the fluid. Initially we expect a hollow-sphere value of 
roughly I in « (2/3)MR 2 « 600, where we used the mass of the outer shell M = 100. 
A more accurate value is obtained by direct summation over the contributing beads, 
yielding I in = 546. In the late-time regime the moment of inertia is expected to be 
the compact-sphere value 7/j n ~ (2/5)MR 2 « 770, where the total mass M ~ 214 
has been used. This results in relaxation times r rin = 0.34 and r T j in = 0.44. These 
values are roughly consistent with a fit to the data in the interval 0.1 < t < 1.0, which 
yields a somewhat larger relaxation time r r = 0.68 (see Fig. EJ). Given the general 
inaccuracy of the estimates, and the crossover of the Debye relaxation function both to 
short-time ballistic behavior, and long-time hydrodynamic behavior, this deviation is 
not too surprising. 

Similar to the translational motion, the exponential decay is then followed by a 
power-law long-time tail. Theory predicts u(i)/u(0) = (irl / p)(47rz/t)~ 5 / 2 [29 . As we 
have seen before, the long-time tail in the translational case is governed by the short-time 
mass as a result of normalization. Similarly, the rotational tail must be controlled by the 
short-time moment of inertia Jj n , which also governs the mean square fluctuations of the 
angular velocity via the equipartition theorem (u 2 ) = 3/cbT '/ Tj n . For this reason, we can 
determine Jj„ by fitting a t~ 5 / 2 law to the data, which is much more accurate than our 
rough geometrical estimate. The fit results in lj n = 533, which is reasonably consistent. 
If we insert this value into the relaxation time expression, we obtain r r = 0.64, in quite 
good agreement with the data. It hence seems that fluid dragging effects are not yet 
very important for the initial Debye relaxation. 

At longer times, one can again notice a significant finite-size effect: the curves 
obtained in the smaller simulation box depart from the asymptotic power-law line earlier, 
i. e. the long-time diffusion is hindered at the small system sizes. In fact, the power-law 
regime is inaccessible at L = 10 while for L = 100 it extends up to t = 200. This value 
characterizes the interval after which the particle starts feeling its own periodic images. 
The rotational diffusion constant D R is given by the Green-Kubo integral [3U] 



We can evaluate this by again making use of linear response theory: The correlation 
function (uj(t) • 3(0)) is identical to the relaxation function presented in Fig. [5J 
multiplied with the initial value (uj 2 ), which we know from the equipartition theorem 
(see above). Using this approach, we have calculated D R for different box sizes L. In an 
infinite fluid, D R has the Stokes-Einstein-Debye value D R = k B T/(8irr]R 3 ) = 0.58 xl0~ 3 , 
towards which the results indeed converge for L — > oo. Again a 1/L behavior (but 
weaker than for translational diffusion) is observed, as shown in Fig. 0] It should be 




(5) 
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Figure 6. Center of mass velocity autocorrelation functions for a neutral and a 
charged colloidal sphere. The Green-Kubo integral function is also shown as solid 
curve for Z = and dashed curve for Z = 10. 

noted that the accuracy of the data is slightly hampered by the fact that we had to 
use a somewhat arbitrary criterion for cutting off the integral function, which does not 
arrive at a constant value but rather ends up with a linear increase originating from a 
small nonzero constant angular velocity in the long-time limit (as a manifestation of the 
finite-size effect). 

Finally, we have started to study the effect of charge on the self-diffusion of the 
colloidal particle. We performed simulations of a sphere with central charge Zq = 10 
in an LB box of size L = 40. Ten counterions of charge —1 were also added. The 
Bjerrum length was set to 2. Technically, the simulation ran without any problems 
just as well as for the neutral system. In Fig. |U]we compare the decay of the velocity 
autocorrelation function for the neutral and charged spheres. The difference between 
the two is not detectable within our error bars at this charge. We cannot say yet much 
about the influence of the charges on the dynamic properties; we do however expect that 
for strongly charged systems the Debye layer will effectively increase the hydrodynamic 
radius and slow down the diffusion. More work needs to be done to resolve this issue. 

4. Summary 

In summary, we introduced and tested a new model for simulating colloidal dynamics 
with inclusion of the hydrodynamic effects on the level of the lattice-Boltzmann 
equation. The suggested "raspberry" -like colloidal object exhibits the essential features 
of the diffusion of a spherical colloidal particle. For the first time we combined a model 
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containing explicit charges with accurate treatment of the hydrodynamics. We expect 
this method to have some advantages over the previously applied scheme [T7| which 
builds upon the model of Ladd |Hj and has the problem of charge leakage. In our method: 
(i) the stick boundary conditions are realized not strictly, i. e. the particle surface is 
softly bound to the fluid; (ii) the method is easily adaptable to the multiple timestep 
method where the ionic action on the colloid can be averaged out somewhat within the 
integration step for the solvent. We should note that such a hybrid simulation requires 
still some computational effort — two million hybrid steps for the colloid consisting of 
100 beads in a box of 20 x 20 x 20 LB cells take about 12 hours on a 2 GHz Pentium 4 
processor, while it is expected to run a few times faster on more sophisticated hardware. 
Here, one step is a MD step of the total system of beads, plus one LB update of the whole 
lattice, where the latter part completely dominates the CPU effort. Some improvement 
is possible by not updating the LB fluid at every MD step, as it was done here. For 
dense systems, where the MD part is no longer negligible, one should also optimize the 
interaction potentials. Finally, we would like to mention that the algorithm can be fully 
parallelized and we are working on the parallel version. 
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